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Abstract 

Lees-Edwards boundary conditions (LEbc) for Molecular Dynamics 
simulations Q are an extension of the well known periodic boundary con- 
ditions and allow the simulation of bulk systems in a simple shear flow. 
We show how the idea of LEbc can be implemented in lattice Boltzmann 
simulations and how LEbc can be used to overcome the problem of a max- 
imum shear rate that is limited to less then 1/L y (with L y the transverse 
system size) in traditional lattice Boltzmann implementations of shear 
flow. The only previous Lattice Boltzmann implementation of LEbc Q 
requires a specific fourth order equilibrium distribution. In this paper we 
show how LEbc can be implemented with the usual quadratic equilibrium 
distributions. 



1 Introduction 

In 1972 Lees and Edwards published a seminal paper (IJ describing an extension 
of the periodic boundary conditions for Molecular Dynamics that allows the 
simulation of a bulk system in a simple shear flow. The uniform shear flow 
steady state reached is computationally convenient because it reduces finite size 
effects. If moving solid walls are instead introduced to produce a shear flow, 
the spatial inhomogeneities induced close to the wall limit the spatial region 
within the simulated system that can be used to study the bulk behavior of a 
system under shear. LEbc ensure that the system is spatially homogeneous, so 
that bulk-like behavior is recovered on smaller simulated systems. For lattice 
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Boltzmann methods there is the additional problem that the maximum shear 
rate is limited by the existence of a maximum shear rate of order 1 / L y with L y 
the transverse system size. This is required to ensure that all fluid velocities 
remain small compared to unity in lattice units S. 

The original LEbc were implemented for a system of particles, with positions 
{ r xir y ,r z ) and velocities {v x , v y , v z ) within a box of dimension (L x , L yi L z ). Lees 
and Edwards showed that by implementing a new kind of boundary condition in 
a MD simulation a state with a uniform shear velocity profile u = jyex (where 
e x is the unity vector in x-direction) could be achieved. The method consists of 
a simple recipe that imposes periodic boundary conditions on particles leaving 
the simulation box in the directions perpendicular to the velocity gradient e y . 
Particles leaving the box in the e y direction (r y > L y ) will be reinserted with 
their x-velocity increased by the shear velocity Au = jL y e x added at a position 
given by their periodic image displaced by the time-dependent offset d x — tAu x , 
t being the time elapsed from an appropriate origin of times. If the particles 
leave the box in the —e y direction (r y < 0) they reappear at the position of their 
periodic image displaced by —d x and have the shear velocity Au subtracted from 
their x-velocity. 

We can describe the boundary conditions by providing the new positions 
r' and velocities v' of the particles after the particles have been moved by the 
dynamics of the algorithm 

(r x mod L x ) + d x r z > L z 

(r x mod L x ) <r z < L z 

(r x mod L x ) — d x r z < 

(r y mod L y ) (1) 




(2) 



This is the best known way of implementing boundary conditions that allow the 
simulation of a bulk material under shear. 

This method works well for all particle-based simulation methods. The aim 
of this paper is to adapt this method for lattice Boltzmann models, where the 
fundamental quantity is a distribution function, rather than a set of particles, 
which is defined on a unit lattice. In this case, the generalization of LEbc is not 
straightforward since the displacement d x will in general not correspond to a full 
integer multiple of the lattice-spacing and it is not clear how to implement the 
required velocity shift in a method with a fixed set of velocities. In an earlier 
paper [0] one of the authors introduced an interpolation method to deal with 
the first problem and a partial Galilean transformation which was used to add 
momentum. This approach required the use of an extended equilibrium distri- 



2 



bution function. Hence, it was not applicable to standard lattice-Boltzmann 
models. 

In this paper we develop an alternative approach to impose Lees-Edwards 
boundary conditions in standard lattice-Boltzmann models. It is based on ap- 
plying a Galilean-transformation and imposing the appropriate change of mo- 
mentum across the moving plane. Additionally, we must displace the densities 
appropriately with an interpolation scheme. The derived method of Galilean 
transformations is quite general and we we can re-derive the correct boundary 
conditions for a moving solid wall in contact with the fluid which were originally 
introduced by A. Ladd Qj. 

After describing the lattice-Boltzmann model, we will discuss in detail the 
implementation of Lees-Edwards boundary conditions, and we will subsequently 
validate the method by examining implementations for binary fluid mixtures in 
two dimensions. 

2 Lattice Boltzmann model 

The lattice Boltzmann method can be viewed as a simple discretization of the 
Boltzmann equation where space is represented by a lattice and the one-particle 
distribution functions (which we will refer to as 'densities') /i(x, t) exist at 
each node x of the underlying lattice, for only a small number of velocities, 
{vi}, that (with unit time step) usually connect a small subset of lattice points. 
These densities evolve following a discrete-time dynamics. In each iteration, 
the densities are moved along their corresponding velocity vectors during the 
streaming step to the appropriate neighboring site, while in the subsequent 
collision step the set of densities at a given node are relaxed toward equilibrium. 
This relaxation step is performed ensuring mass and momentum conservation 
at each node. The energy is usually not conserved (although energy conserving 
schemes also exist ||) but it is possible to introduce a nominal temperature, that 
can be viewed as the effect of an underlying thermostat. This temperature is 
introduced as a parameter of the equilibrium distribution function (see below). 
The evolution of the densities in each time step can be written as 

/i(x + v< , t + 1) = /i(x, t) + ^ (/P(p(x, t), u(x, t)) - /<(x, *)) (3) 

where r is a relaxation parameter, which in the hydrodynamic limit of this 
equation is related to the viscosity of the fluid. Eq.(||) can be viewed as the 
discretized version of BGK approximation of the Boltzmann equation. The 
linearized collision operator can be generalized to allow for several relaxation 
parameters (and hence the existence of different viscosities for the fluid) In 
numerical applications the above evolution equation is usually split into two 
subroutines, a collision and a streaming step. The collision step is local 

/i(x,t) = /((x,t) + - (/.«V(x,i),u(x,i)) - ffat)) (4) 
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and the streaming step moves densities along their associated velocity vectors 



f i (x. + v i ,t + l) = f i (x,t). (5) 

In the above equations, is the equilibrium distribution function that deter- 
mines the thermodynamics of the fluid, while p and pu are the local mass and 
momentum densities. These are the macroscopic quantities of interest, and can 
be obtained from the densities fi as appropriate moments 

i i 

The choice of the equilibrium densities is one of the key ingredients of the model. 
They are chosen to recover the equilibrium density and momentum 



P 



=E/° /™=£/N (7) 



and the appropriate pressure tensor in equilibrium 

Va/3 = £ -ft( Via ~ U a){ v iP ~ u fi) ( 8 ) 
i 

where V is the sum of the thermodynamic pressure tensor P and an additional 
term that is required for Galilean invariance of the system j(| : 

r — 0.5 

Vap = P a (3 H ^— {u a dpp + u fj d a p) (9) 

Using the thermodynamic pressure tensor P is one of the standard ways to 
impose the desired thermodynamic behavior of the lattice-Boltzmann fluid |Q. 
In the simplest case, the pressure tensor for the ideal gas is P a p = p^8 a p. We 
need an additional constraint for the third moment 

^ ftV la Vij3V l7 = ^(Mq(5/3 7 + Up5 aJ + UySap) (10) 

i 

to ensure that the anisotropy of the underlying lattice does not affect the effec- 
tive behavior of the model in the hydrodynamic limit. It is possible to fulfill 
these constraints by introducing equilibrium densities that depend quadratically 
on the lattice velocities, {vi}. Depending on the geometry of the lattice, there 
is still some freedom in the choice of the equilibrium densities. 

A Taylor expansion to second order in the derivatives [[| of the zeroth and 
first moments in v of the lattice Boltzmann equation (|^) gives the macroscopic 
equations which express mass and momentum conservation. From (^) we have 

£- T (3t+V i Vr/ i = -(/?-/;) (11) 
n=l 
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and wc can approximate 



Si = /° + r(d t f? + v ia d a f?) + 0(d 2 ). (12) 

The zeroth velocity moment of (^) up to 0{d 2 ) gives the continuity equation 

d t p + d a (pu a ) = (13) 

and the first velocity moment gives the Navier Stokes equations (with <9 a uuu = 
0(d 3 )) 



pd t u a + pupdpu a = -dpPap + dp 



pv ( d/3u a + d a Uf3 - -^<9 7 u 7 <5 Q/3 



(14) 



where the viscosity is v — (t — 0.5)/3 and D is the number of spatial dimensions. 
In order to simulate a mixture of two components A and B we follow OrlandiniQ 
and define a second lattice Boltzmann equation 

gi (x + Vi,t+l) = gi {x,t) + — [^(p(x,t),c£(x,i),u(x,*)) -ft(x,t)] (15) 
where the density difference (concentration) 4> = pa — Pb is defined as 

= (16) 

i 

and is the only additional variable conserved in the dynamics. The total density 
is now p = pA + Pb and the concentration couples back into the momentum 
equation, eq.(|l4j), via the pressure tensor P a $ which is now a function of both 
p and (f>. The relaxation step for the dynamics of gi is dependent on which is 
related to the concentration diffusivity in the hydrodynamic limit. The equilib- 
rium densities gi are determined, again by imposing that the local equilibrium 
concentration is recovered (from eq.(|l6|)). The additional constraints on the g® 
are 

<j>Ug = / J 9jVia, M a p — 5i \ v ia ~ U a )(vip - Up) (17) 
i i 

Analogously to the pressure tensor, M. is given by a correction term for Galilean 
invariance [|| and the chemical potential p: 

M a f) = Tp8 a p + - - ° 5 (uadpcf) + upd a <f>) (18) 

The desired thermodynamics of the fluid is determined by p which must be 
consistent with the chosen P a p. These two quantities are normally derived from 
a given free energy || fjj. In this way one ensures their consistency, and can fix 
beforehand the appropriate thermodynamic behavior of the binary fluid model. 
It is enough to consider equilibrium densities that depend quadratically on 
the velocities to impose the previous requirements. 
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Again, we recover the hydrodynamic limit of eq. jl5|) performing a Taylor ex- 
pansion of to second order in spatial gradients. In this case we get a convection- 
diffusion equation for the concentration 



d t (t> + d a {(/)u a ) = d a < (T0 - 0.5) 



9/3 (Tp(p, 4>)8 a p) + -dpP a/3 



(19) 



This model has been used extensively for the simulation of two-phase flows 
[|[ [|, |(| 0, |).The following derivation of the LEbc presented in the next section 
is, however, even more general and can be applied to most lattice Boltzmann 
methods. 



3 Lees Edwards' boundary conditions 

The implementation of Lees Edwards boundary conditions for lattice Boltzmann 
will follow the same steps as in its original implementation: we will apply pe- 
riodic boundary conditions to the distribution function in the two directions 
perpendicular to the velocity gradient; and in the third, we will assume that 
the periodic system is moving as a whole with a velocity Au with respect to 
the lattice. In a second step, we will have to displace the densities fa (and gi 
for the binary model), taking into account that we will have to interpolate the 
densities from the off-lattice positions to which they stream, onto the regular 
lattice. So this boundary condition will be applied after the collision and before 
the streaming of the densities. 

Let us first focus on the velocity transformation. Since we are dealing now 
with a distribution function, the way to implement the velocity jump when 
crossing the system boundary in the velocity gradient direction will be to per- 
form a Galilean transformation, i.e. if the distribution function in the system 
has momentum pu, its image node will have momentum p(u + Au) if we move 
in the direction of the gradient, and p(u — Au) when we leave the system in 
the direction opposite to the velocity gradient. In this way, we will generate a 
linear velocity profile with shear rate r yAu/L y . In order to get the momentum 
transfer, we need to calculate the difference between a density fa at the system 
velocity and the one after the Galilean transformation has been applied. All the 
macroscopic variables and their derivatives will be the same, since the image 
system is moving as a block at constant velocity. So, using eq.(|l^), where we 
retain only the it-dependence of the /'s, we obtain 

A/, = /,(u + Au)-/,(u) 
« /P(u + Au)-/P(u) 

-r{9 t [/P(u + Au) - /°(-u)] + v ia d a [f?(n + Au) - /°(u)]}(20) 

where Afi is the change in the density moving in the i direction due to the 
Galilean transformation. Note that the transformation will only be applied to 
the densities crossing the LE boundary but not to the rest of the densities at 
the same node. 



G 



The terms linear in r are of order O(d) and their zeroth and first moments 
are 0(d 3 ). These terms are therefore small and enters only into the pressure 
moments. Therefore we can neglect them. We obtain for the boundary condition 

/j = / i + /°(u + Au)-/°(u). (21) 

where f[ means the density after the application of the boundary condition. 
We now need to specify the equilibrium distribution. We will assume the usual 
quadratic distribution function which can fulfill the conditions @) and (§): 

f° = pa l + pa\v ia u a + a\v ia VipTl a p + a3tr(II) (22) 

where Hap = P^a + P u aUj3 and we obtain 

/?(u + Au)-/°(u) = p[a\v ia Au a 

+a2V ia Vip(u a Aup + UjjAu a + AllaAltp) 

+4(2u.Au+ |Au| 2 )] (23) 

If we have a binary mixture, the same transformation will apply to the densities 
gi related to the concentration. The equilibrium distribution g® can be written 
in the same form as eq.(p^) if we replace p by <fi and P a p by Tp,5 a p. In the 
previous expression T is a factor which, together with the relaxation parameter 
T0, will determine the diffusivity of the concentration; 5 a p is the Kronecker delta 
function. 

Now that we have defined a suitable Galilean transform, we have to split the 
densities into an untransformed part and a transformed part that will stream 
over the boundary, i.e. the densities at the top and bottom boundaries with 
associated velocities Vi y > and Vi y < respectively if we take e y as the 
direction of the shear gradient. It is important to ensure that the two subsets of 
velocities conserve the node density after the Galilean transformation has been 
applied. This requires 



p = E /<+ E u 

Viy>0 Viy<0 

= ]T [h + /°(u + Au) - /°(u)] + E fi 

V iy >0 V iy <0 

= £ [/°(u + Au)-/°(u)]+E/* 



V iy >0 

^ = ^ (4 + ^) [ / o (u + Au) _ / o (u)] 
1 



2 (Pyy + P u y u y ~ P yy ~ pu y u y ) (24) 

where Y] v . >(<)o i s performed over the velocity directions with a positive (neg- 
ative) component crossing the shear boundary. Eq. (E4J) is fulfilled if we require 
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that the velocity subset that characterizes the lattice Boltzmann model involves 
only velocities that do not span beyond the closest site layer, i.e. the factor 
i v i y + u iy)/2 is 1 for velocities with v y > and otherwise. This is the case 
for most standard velocity sets and for the velocity set we will consider here 
(see eq.(^)) in particular. To derive this result we have also used that the 
velocity-shift is parallel to the boundary, i.e. Au y = 0. 

It is worth remembering that the total momentum in the direction of the 
velocity gradient in the system must be zero; otherwise there will be a net 
mass transfer across the LE boundaries which would result in an increase (or 
decrease) of momentum in the direction of the velocity gradient. The mean 
velocity will then steadily increase until it exceeds the maxium velocity and 
the simulation becomes numerically unstable. This is a general requirement of 
any Lees-Edwards boundary condition, because the presence of LE boundaries 
breaks the momentum conservation unless the mass transfer is exactly balanced 
in both directions. 

The idea of the Galilean transformation that we have derived also has more 
general applications and can be used to derive the boundary conditions for a 
solid moving wall. At a solid wall, the velocity of the fluid must be that of 
the wall. This no-slip boundary condition for a stationary wall is implemented 
through the 'bounce-back' method In it, one replaces the streaming step of 
eqn.(||) for the densities that would have been streamed through the wall with 

/ i (x,t+l) = /_ i (x,t) (25) 

where V-i = — M,-. To generalize this bounce back mechanism for a moving 
wall we first perform a Galilean transformation of the densities into a reference 
frame in which the wall is at rest. We then perform the bounce back and then 
transform the densities back with an inverse Galilean transformation. We then 
obtain for the bounce back off a wall moving with velocity U 

/i(x, t + 1) = /_i(x, t) + 2pa\v ia U a (26) 

which is a result first derived with a different method by Ladd^, eq.(3.3)]. 
For this particular case, the change in the distribution function is linear in the 
velocity jump, rather than quadratic, as in the general case (cf. eq.|23|) 

To complete the implementation of Lees-Edwards boundary conditions we 
also need to move the densities by the displacement of the image system, d x — 
tAu x , which is not going to be an integer in general. We must then split the 
displacement in an integer part d\. and a real part with < d x < 1. We 
then shift the densities that will stream over the boundary with Vi y > and 

1 According to the discussion of the previous paragraph, the density at both sides of the 
boundary is not modified by the bounce in the case of a moving parallel solid wall that slides 
along one of the boundaries of the system. For more general situations, this is not the case. 
For example, in the case of a moving rigid object embedded in the fluid, its velocity is not 
parallel to the boundary. In this case, the density inside the rigid object will be modified as 
a function of time. However, the presence of a rigid interface separating the inner from the 
outer motion of the fluid makes it possible to correct for these deviations [TL if] . 
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Table 1: The coefficients for the D2Q9 model. There are two free parameters 
that we set to k = I = 1 for convenience. 

Vi y < respectively by 

/'(x± t) = (1 - d«)/(x + 4, t) + d*/(x + 4 ± 1, t) (27) 

where are the positions of the lattice points directly below or above the 
LE boundary and this shift and interpolation is applied to densities /j with 
Viy = ±1. That is, we have performed a linear partitioning of the moving 
density between the two corresponding nearest nodes. If necessary, it might be 
possible to improve upon this linear interpolation. This simple choice, however, 
already gives very accurate results. 

Lees Edwards boundary conditions for lattice Boltzmann were first derived 
by Wagner & Yeomans Q whose approach also required 

A/ i = /P(-u)-/°(u) (28) 

but the authors also insisted on a heuristic transformation rule where the amount 
of velocity added to the densities that stream over the boundary should corre- 
spond to Au. In particular they first calculated a separate x-velocity for the 
particles that crossed the LEbc and then they required that this x-velocity 
should be independent of the y-velocity in equilibrium. In order to fulfill this 
additional requirement they needed a more complicated equilibrium distribution 
(a polynomial fourth order in u) . While this additional requirement is perfectly 
sensible in the limit of the continuous Boltzmann equation our results show that 
it is not always helpful to consider the densities to be those of real particles and 
that it is perfectly feasible to have a Lees Edwards boundary condition for a 
lattice Boltzmann method with a traditional quadratic equilibrium distribution. 

4 Numerical Implementation 

We implemented a lattice with dimensions L x ,L y and the boundaries of the 
simulation lattice have periodic boundary conditions. We implemented a D2Q9 
model, a two dimensional model with velocities 

-{ffl.Q.(-j).(?).(-!).©-(-;)-(=;)-(-;)} « 

and with the coefficients given by Table 1. 
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4.1 Multiple LE planes 



One of the limitations of Lattice-Boltzmann models to study fluids under shear 
arises because the maximum velocity cannot be larger than one lattice spacing 
per time step (and in practice must be ~ 0.1). This gives an upper bound of 
order X/L y for the shear rates that can be reached with the method described 
so far (one cannot simply increase the lattice size while keeping the shear rate 
constant). This suggests that large systems with Lees-Edwards boundary con- 
ditions should be restricted to small shear-rates. We can now overcome this 
limitation using the present method by introducing a larger number of Lees- 
Edwards shear planes as one increases the L y dimension of the simulation. We 
will still require, of course, that the velocity in each of the sub-lattices separated 
by the LEbc are no larger than about 0.1. 

We will allow for N LE Lees-Edwards shear planes at almost equidistant 
positions 



where [x]j is the largest integer smaller than x. The system is then effectively 
cut into N LE — 1 sub-lattices that will move relative to each other with the 
associated jump- velocities U k . The average shear rate for the system is therefore 



and the stationary solution due to this set of LE shear planes is a linear shear 
profile of the form 



where j/q is the position of zero velocity for each sub-lattice. 

At this point it may appear that the introduction of the shear-planes does 
introduce a local velocity profile and biases the flow but, in the absence of 
errors, the only the global shear rate 7, defined in eq.(|3l|) matters to the flow. 
This is easily seen if one considers a Molecular Dynamics simulation where 
the introduction of several shear planes is equivalent to describing the different 
sublattices in different frames of references, with no effect on the real dynamics. 
Because of this there is, of course, no merit in introducing more than one LE- 
plane (which establishes 7) in MD simulations. The only limit for the use of 
many LE-planes in lattice Boltzmann lies in the errors induced by the additional 
internal boundary condition and we discuss them in the following sections. 

4.2 Transient start of shear and stationary profile 

As a first test for the new boundary conditions we will examine the accuracy 
with which the algorithm simulates the start-up of a shear in a single component 
fluid. We consider a system with one LEbc and the initial velocity is zero. This 



y = Y k =[(k-l/2)L y /N LE ] I forfc = {l,...,iV iB } 



(30) 



7 = 




(31) 




(32) 
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Figure 1: Start-up of a shear profile. Solid lines are given by eq. ( |36| ) and 
the symbols represent simulation data. The inset shows the difference of the 
theoretical and simulated velocities. 
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corresponds to a system in which there is a periodic stack of slabs of fluid of 
width L y each of which moves with a different velocity. The velocity of the nth 
slab is given by u x = uAu. The system will then exhibit viscous friction and 
the profile will converge toward a linear shear-profile. The density and pressure 
in this system remain constant and the velocity only has an x-component and 
its value only changes in the y-direction. The Navier-Stokes equation for this 
system gives 

d t u x = dyU x (33) 

i.e. a simple diffusion equation for the momentum. The solution for this equa- 
tion for L y — oo is well know |L2], §52] and given by 

u x (y, t ) = erf (^=) (34) 

where 

erf(x) = 4= / ( 35 ) 
V n Jo 

Obtaining the solution for finite L y is simply a linear superposition 

°° k 

«*(!/,*) = J2 W\ 

fc=-oo 1 1 

which can be conveniently evaluated numerically. The agreement between the 
theoretical prediction and the measured velocity is very good as can be seen in 
Figure H. The error shown in the inset of Figure [j] is largest for the early time 
simulations but is always below 1% of the maximum velocity. 



1-erf 



kL, 



2Jvi 



(36) 



4.3 Shear in a two component model 

The previous simulation was effectively a one-component simulation and it did 
not test the effectiveness of the interpolation scheme. Unfortunately, it is not 
easy to find two-component flows under shear for which analytical solutions 
are known. There is, however, a method to test the errors introduced by the 
LEbc for a two-component system very accurately. We can choose to have two 
equidistant LEbc with equal but opposite velocities so that the effective shear- 
rate is zero. The resulting system is then equivalent to a system at rest. We will 
consider two periodic stripes of A-rich and B-rich material at rest and simulate 
it with two LEbc with opposite velocities. 

In figure || we show the results of a simulation on a 40x60 lattice with a 
stripe that was 20 lattice spacings wide. There were two LEbc with velocities 
of 0.005 and -0.005. The figure shows a section of 20x20 after 20000 iterations 
(d x = 0) where the stripe interface overlaps the LE boundary. The full system, 
of course, comprises of two interfaces in the x-direction and two LEbc in the 
y-direction and there the same error-patterns appear by symmetry. 

Without errors this system would show perfectly straight contour lines and 
velocities of 0.005 in the upper half of figure |^(a) and -0.005 in the lower half. 
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(a) (b) 

Figure 2: Steady state profile of the stripe morphology. In (a) the real velocities 
are shown and the maximum velocity corresponds to 0.0051. In (b) the deviation 
from the theoretical profile is shown and the maximum velocity here is 0.00013. 
Eight contour lines are also shown indicating the width of the interface. (Sec 
text for details) 



Since the velocities are very nearly that, we have to examine the difference of the 
real velocity and the predicted velocity shown in figure ||(b) . We see that this 
error corresponds to a deviatoric velocity profile that flows into the interface at 
the LEbc. We attribute this error to the fact that we use the linear interpolation 
from eq. (p7|) . The effect of this interpolation is to smear out the interface and 
make it wider than the equilibrium profile. There will then be a Marangoni flow 
that tries to relax the profile to its equilibrium shape which we assume is the 
major contribution to this velocity profile. A close examination of the interface 
reveals a slight widening around the LEbc. The maximum amplitude of the 
error in the velocity is 3%. In future work we plan to examine if we can reduce 
this error by using an interpolation method that will conserve not only the mass 
but also the derivatives and prevents the additional diffusion introduced by the 
LEbc. Nonetheless, the accuracy found above is comparable to other sources of 
error (e.g. anisotropy) present in lattice Boltzmann as currently used. 

As a final benchmark of the performance of LEbc's, we have studied the 
shape of a drop in a uniform shear flow generated by a single LE plane. In 
order to check for artifacts induced by the presence of the boundary, we have 
compared a simulation where the drop does not cross the boundary, and a 
situation where the LE plane bisects the drop. In Fig.|| we show the comparison 
for two different capillary numbers Ca -defined as the ratio between viscous and 
interfacial forces, Ca — vRrf/a, with R the radius of the drop and a its surface 
tension. The interface is located in the simulation through an interpolation of 
the corresponding nearby nodes where the order parameter changes sign, which 
introduces a small source of error. The drop that is bisected by the boundary 
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(b) 



Figure 3: Steady shapes of a drop in a Couette flow for different capillary 
numbers. In the open circle situation the drop does not intersect the LE plane. 
The filled symbols correspond to the simulation where the LE plane bisects the 
drop, a) Ca = 0.012, b/ Ca = 0.5. The interfacial width in both cases is of the 
order of 3 in lattice units. 



will have a larger averaged velocity. Since we are comparing the drop shapes 
at equal time, we have displaced the drop accordingly to make them overlay. 
In the figure one can clearly see that the intersection of the boundary does not 
affect the shape of the drop, within the numerical accuracy already mentioned 
in the previous paragraph. 



5 Conclusions 

In this paper we have described how to implement Lees-Edwards periodic bound- 
ary conditions for lattice-Boltzmann models. Lees Edwards boundary conditions 
are a useful for the simulation of bulk systems under shear without the compli- 
cating influence of walls. In this way, finite size effects are diminished, which 
has converted the uniform shear flow in an attractive steady configuration to 
analyze computationally. 

The implementation we have described follows the basic steps of the origi- 
nal paper and is compatible with the standard LB models, that assume an 
equilibrium distribution for the densities which is quadratic in the velocity. The 
numerical examples described show that the method is quantitatively accurate 
for single and even for two component systems. 

We have also described how to avoid the limitations on the attainable shear 
rates that are imposed by the underlying lattice. By adding a number of shearing 
planes, steady states of large shear rates can then me reached (subject only to 
the intrinsic stability limitations of any lattice-Boltzmann schemes). Because 
of the errors introduced by the interpolation scheme, it is not advisable to 
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introduce more LEbc's than necessary. One could think of a system where 
the whole lattice is subject to LEbc but this would now lead to a non-isotropic 
diffusion. We are planing to improve the interpolation scheme in order to reduce 
the artificial diffusion. 
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